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We study localized plasmons at the nanoscale (nano-plasmons) in graphene. The collective 
excitations of induced charge density modulations in graphene are drastically changed in 
the vicinity of a single impurity compared to graphene's bulk behavior. The dispersion of 
nano-plasmons depends on the number of electrons and the sign, strength and size of the 
impurity potential. Due to this rich parameter space the calculated dispersions are intrinsically 
multidimensional requiring an advanced visualization tool for their efficient analysis, which 
can be achieved with parallel rendering. To overcome the problem of analyzing thousands 
of very complex spatial patterns of nano-plasmonic modes, we take a combined visual and 
quantitative approach to investigate the excitations on the two-dimensional graphene lattice. 
Our visual and quantitative analysis shows that impurities trigger the formation of localized 
plasmonic excitations of various symmetries. We visually identify dipolar, quadrupolar and 
radial modes, and quantify the spatial distributions of induced charges. 
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1. Introduction 

Graphene is a single-layer atomically thin allotrope form of carbon, where carbon 
atoms are arranged in a honeycomb lattice. Owing to its single-atom thickness 
and unique electronic properties it is a promising material for technological use [1] . 
Although graphene was originally thought to be a pure system, recent scanning 
tunneling microscopy work shows that it is a disordered system [2J. It has lattice 
defects, vacancies, ripples, dislocations, magnetic impurities, etc. Particularly there 
are charge puddles caused by chemical adsorption or molecules trapped between 
the graphene lattice and substrate [3] . Previous studies of impurities have focused 
on understanding their effects on the ground state electronic properties of graphene 
[2t However, an equally important aspect is the study of the consequences 

of these charged impurities on the dielectric response for future electronic devices. 

Several studies explored plasmonic excitations in pristine graphene fTSHTSj . So far 
almost no work focused on localized plasmons (nano-plasmons) in graphene. These 
excitations have the ability to concentrate an electric field (visible photons) at the 
nanoscale regime, thereby beating the diffraction limit. The potential technological 
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applications of this effect are vast. For example, it is possible to improve solar 
cell efficiency by increasing the cross section of the absorbing material without 
requiring a thick layer of it [19j. Another potential application is the new sub- 
field of photonics called plasmonics, where plasmons substitute photons in the 
information transport extending photonics into the nanoscale world [20j . In biology 
the potential applications include labeling and imaging, optical sensing and photo- 
therapy treatment of cancer [21j. 

Here we study nano-plasmons, that is, localized plasmonic excitations in pristine 
and disordered graphene. When plasmons are localized the translational invari- 
ance of the lattice is lost, hence analytical calculations become impractical. A 
self-consistent quantum-mechanical approach to nano-plasmons in graphene in the 
presence of impurities has been developed, which accounts for the non-locality of 
the dielectric response function (SUES]. Here we use a generalized version of this 
approach, where the polarization operator is diagonalized, providing information 
about individual plasmonic excitations, including spatial profiles and local spectral 
density of states. The calculation is performed for a single impurity at the center 
of the carbon hexagon. We study the impurity effect as a function of the sign, 
strength, and size of the impurity potential, and the doping level of graphene. The 
interactive use of three-dimensional (3D) visualization allows us to identify readily 
which plasmonic excitations are localized vs. delocalized. Our visualization of the 
poles of self-consistently calculated polarization operators shows the existence of 
nano-plasmons around naturally occurring impurities in graphene, which confirms 
that graphene is an intrinsic plasmonic material [24j. 

The paper is organized as follows: In Section 2, we introduce the lattice model 
used to calculate the plasmonic excitations in graphene. In Section 3, we describe 
the needs for an advanced 3D visualization tool to analyze complex multidimen- 
sional data sets and discuss the computational requirements. The results of the 
plasmonic modes are presented in Section 4. In Section 5, we perform a quantita- 
tive analysis of the symmetries of the various nano-plasmonic modes and conclude 
in Section 6. 



2. Lattice model 

The electronic structure of graphene is described by the Hamiltonian H — + V ^ 
where is the kinetic and potential energy, while V is the Coulomb interaction 
between electrons. The interaction is treated as a perturbation and will be consid- 
ered later in the Random Phase Approximation (RPA). The single-particle part of 
the Hamiltonian on a lattice is 

<ij> i i 

where < ij > denotes nearest neighbor sites; cj and q are operators that create and 
annihilate electrons at site z; t = 2.7 eV is the tight-binding hopping parameter, and 

fi is the chemical potential (/x = for undoped graphene). Ui = Uq exp (^ ~^^2a^^^ ) 
is the change in the on-site potential for atom at site due to the presence of an 
impurity at xq. Uq is the strength of the impurity potential and a determines its 
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spatial spread. The matrix structure of the non-interacting Hamiltonian is 



'■■ Ui-i + n -t 
-t Ui + fj, 







-t Ui+i + iJ, 



(2) 



The diagonalization of the non-interacting single-particle Hamiltonian gives 
eigenstates |^^^) and eigenvalues which satisfy J2j I'i^^)^ = |^^^)^. 
Each state has an occupation number given by the Fermi function = 
[exp(-^) + l]-i. 

In order to compute the polarization operator due to an induced charge, we 
define the tensorial matrices An and AE in the basis of the eigenstates 



(3) 



Using the superindices I = a® ^ and J = ^ ® 5^ which reshape a matrix into a 
vector in a column- wise fashion so while a, /3 G 1, . . . , n one has / G 1, . . . , n^, the 
above tensors can be written as diagonal matrices 



An 



nx — n 



nx - m 



(4) 



and 



AE = 



E^^-E^^ 



— '^^f+l. 



(5) 



The polarization operator is a density-density response function. Summing the 
product of non-interacting Green's functions over Matsubara frequencies Um one 
obtains the non-interacting polarization operator 



= s, 



(6) 



Using the superindices / and J, one can write as the product of diagonal 
matrices 



{u;) = An{AE-u;I)-\ 



(7) 



which again is a diagonal matrix. 

Next we include the Coulomb interaction V for electrons on the lattice. The 
interaction Hamiltonian is written in terms of creation and annihilation operators 
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in the lattice basis as 



ijlm 



where the only non-zero components are 



We can write V in the basis of the eigenstates, where it will be a full tensorial 
matrix. Therefore V can also be represented as a matrix in the superindex space 

Within the RPA calculation the polarization operator of the interacting system 
is given by 

u{oj) = {u) (1 - vu^ {u)y^ . (10) 

Using Eq. Q we can write 

n (cj) = An {AE - VAn - uiy^ . (11) 

The nano-plasmonic excitations correspond to the poles of the polarization opera- 
tor. Hence the modes correspond to charge density oscillations with amplitude 

5p/ = 5^An/j(/)j, (12) 
J 

such that 

{AE -ul- VAn) = 0. (13) 

The matrix AE — AnV is diagonalized giving the eigenvalues uo"^ and eigenvectors 
satisfying ^j{AE — AnV)u(pj = cj^(f)j- The polarization has poles at = cj^ 
for each eigenvalue ou^ of AE — AnV. One can then find the spatial profile of the 
mode with energy ou^ through 6pj = J2j ^^u^j- going back to the lattice 
basis and neglecting the orbital overlap we get the local charge oscillation 

SPI = 5Pa/? ^ Spfj ~ 5pi (14) 

The amplitude of the induced charge density oscillation at site i is the vector 
6pf = Spf-. Of course the induced charge will oscillate with frequency cj^, that is, 
6p^{t) = 6p^exp{-icjH). 



We use a standard numerically exact diagonalization method to solve Eq. (13). 
This allows us to find all the poles of the polarization operator and, thus, we obtain 
the complete information of the plasmonic excitations on the graphene lattice. For 
a lattice having N sites, the matrix size of 11 becomes = N'^ x N'^. The number 
of eigenmodes scales as = N/2 x N/2 for N/2 electron-hole pairs at half- 

filling. Therefore the number of data points where the plasmon dispersion needs to 
be visualized scales as which becomes rapidly a big number even for small 

lattice sizes. 
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3. The need for advanced 3D visualization 

Visualization plays an important role in scientific data analysis with a growing 
need for being able to visually analyze complicated data sets in three dimensions 
(in real space) and more (in parameter space). The analysis of such data is always 
complemented by the ability to effectively interact with the data. Interacting with 
multidimensional data sets in a 3D environment (stereo display) adds an additional 
advantage to navigating and seeing through complex data sets. Although this ad- 
vanced visualization feature is not always necessary for the data analysis. The need 
for visual analysis of large data sets (larger than the size that a typical computer 
with a graphics card can process) is growing. The problem is commonly tackled 
by using distributed parallel processing on multiprocessor units. Here we discuss a 
specific example of analyzing a multidimensional, medium-sized data set (on the 
order of Megabytes). The data set has 200,000 to 300,000 points in three dimen- 
sional space. In this paper, we show the efficiency of a small-scale parallel rendering 
computer cluster as a possible solution. The data is obtained from a calculation of 
the spatial distribution of plasmonic excitations (induced charges on lattice sites) 
in pristine and impure graphene. Since visual analysis of data combined with quan- 
titative analysis provides a better understanding of the data studied, we present 
specific examples for the quantification of visual information using the open-source 
visuahzation tool ParaView [25j. 

We analyze the dispersion of the plasmonic excitation as a function of energy 
(jj^ ^ uj. Since the graphene lattice is two dimensional, the analysis of the induced 
charge distribution on individual lattice sites as a function of energy uj adds one 
more dimension to the visualization problem. For a given impurity potential, a 
single data set is already four dimensional (4D) in the hyper-space of variables 
(X, Y, cj, 5p(X,Y,Cl;)). The X-Y dimension is fixed by the number of lattice sites 
chosen for the calculation. We fixed it to be = 96 lattice sites. The Z-dimension 
is determined by the number of possible energies uj of the plasmonic excitations. 
For 96 lattice sites the number of unique frequencies oj is about 2300, and the 
number of induced charges as a function of frequency to be visualized is about 
^ 220,000. To specify the induced charge oscillation at each lattice site, 
we need one scalar value and three coordinates. Further, we need at least two 
colors to represent the polarity of the induced charge. In three dimensions it is 
necessary to use a 3D glyph to represent the scalar value (induced charge). We 
use a "sphere glyph" having at least 10 points of resolution for angle and 10 
points for angle [26^ i27j. In the XML unstructured grid data file format [26r.-29j , 
for which the visualization tool ParaView has a reader to input data, the size of 
the raw data file is about 10 Megabytes. After using the 3D sphere glyphs, to 
represent the induced charge distribution, the size of the data set to be processed 
becomes about 1.3 Gigabytes ^ 2 x 10^ x 10 Megabytes. So the memory required 
by the computer to process the data becomes more than 100 times larger than the 
size of the original raw data set. One can use a single-processor computer with a 
high-end graphics card to analyze and post-process a couple of such data sets, but 
many data sets of this size cannot be efficiently processed, let alone attempting a 
real-time or interactive data analysis. Hence using a single-processor computer is 
neither time efficient nor practical and parallel rendering on multiple cores or CPUs 
(central processing units) or multiple GPUs (graphical processing units) becomes 
mandatory. In the remainder of this study we demonstrate how distributed parallel 
rendering on multiple CPUs can solve the problem of large data sets. 

In our case, we study the plasmonic excitations as a function of several parame- 
ters, such as the sign (taking positive and negative sign of Uq)^ strength (magnitude 
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Figure 1. (a) 2D graphene honeycomb lattice with a single impurity (turquoise sphere) at the center of 
one of the hexagonal cells, (b) We show the distribution of induced charges on some of the lattice sites 
of pristine graphene as a function of plasmon energy (third dimension or Z-axis, i.e., direction of skewers 
pointing out of the plane above each lattice site). The size and color of glyphs represent the magnitude and 
polarity of the induced charge, respectively. We show the charge distribution only in the region of interest 
near the impurity site. Here we set the impurity potential Uq = 0. 

of Uq)^ and size (varying a) of the impurity potential, and the doping level (chang- 
ing fi) of the lattice. To understand the interplay of these parameters in determining 
the plasmonic excitations we have to generate many data sets. By varying just one 
of these parameters we introduce an extra dimension to the already 4D visualiza- 
tion problem, namely a fifth dimension. A visualization tool that can handle this 
extra dimension is desirable. For example, in ParaView one can encode the fifth 
variable (or dimension) as virtual time and selectively scan through a data set. 
The selective rendering of the data set, corresponding to different variables, can 
be done from both a graphical user interface (GUI) and script user interface (SUI) 
[261 ET]. The SD-rendered object can be displayed for multiple variable parameters 
allowing a comparative study of the response of the system to these parameters. 
Additionally, a quantitative analysis of the data set can be performed by using spe- 
cial filter functions already provided by the visualization tool or custom developed 
[261 EZ]. 

For time efficient visualization, we used the server-client mode of parallel ren- 
dering in ParaView to analyze and visualize the data sets. We checked the scaling 
of the problem by recording the rendering time for different number of processors. 
We found that it scaled linearly for up to 16 processors (the maximum number we 
tested). The results presented below were obtained by scanning through the data 
sets in the multidimensional parameter space. 



4. Plasmonic modes in graphene 

In Figure [T]^a) we show the geometry of the graphene lattice used for our calcula- 
tions. The black dots specify the positions of carbon lattice sites. The honeycomb 
lattice used in our calculations has = 96 sites. The turquoise colored sphere, 
which sits at the center of a hexagonal cell, represents the position of the sin- 
gle impurity. All results presented here are for a single impurity on a 2D lattice 
with periodic boundary conditions. The properties of the impurity are modeled by 
their effects on the on-site potentials at neighboring carbon sites. Unless otherwise 
specified, the spatial extent of the impurity potential is restricted to the nearest 
neighbor carbon sites. Throughout this work, the energy is measured in units of 
the hopping energy, t = 2.7 eV, between nearest neighbor carbon sites of pristine 
graphene. 
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Figure 2. Distribution of the induced charges, on several lattice sites near the impurity, as a function 
of plasmon energy ujjt. Results are shown for t/o/t = 0.47 (a), 0.79 (b), 1.42 (c), and 2.05 (d), where 
t = 2.7 eV is the hopping parameter, which is used as the unit of energy. The strength Uq of the impurity 
drastically modifies the local response of the lattice. The impurity strongly localizes the induced charges 
near the impurity site. 



4.1. Pristine graphene 

To study the plasmonic excitations in pristine graphene, we set the impurity po- 
tential ?7o = in our calculations. In Fig. we show the distribution of induced 
charges for a small section of the graphene lattice as a function of energy. The 
Z-axis (here and in what follows) represents the energy uj of the plasmonic exci- 
tations. We show the charge distribution only for a small section of the lattice, 
because later on this will be the region of interest for a spatially small-sized im- 
purity potential. The color and size of the glyphs represents the polarity and the 
magnitude of the induced charges, respectively. This visual image will be used as 
a reference image for a comparative study of the effects of a single impurity on 
nano-plasmonic excitations concentrated near an impurity. 



4.2. Effect of positive and negative impurity potential at zero doping 

We study the effect of a single impurity on the spatial distribution of the induced 
charges in the close vicinity of the impurity site for all plasmonic energies. First 
we discuss the effect of an impurity with a repelling (positive) potential. In Fig- 
ure [2] we show results for the impurity potential ?7o/t=0.47, 0.79, 1.42, and 2.05, 
respectively. Comparing Fig. [T]^b) with Fig. [2]^a), one can see that the region clos- 
est to the impurity site shows the largest change in induced charge distribution. 
As the impurity potential is increased the change becomes more pronounced, see 
Figs. [2]^b)-(d), as seen in the increase in the size of the 3D sphere glyphs used to 
represent the magnitude of the induced charges. More precisely, induced charges 
start increasingly to localize on neighboring sites of the impurity. These are the 
sites which are maximally affected by the impurity potential. 

We also study the effect of an attractive impurity potential on the induced charge 
distribution for Uo/t = —0.47,-0.79, —1.42, and —2.05. Similar to the case of a 
positive impurity potential, the most pronounced change occurs on lattice sites 
in close proximity to the impurity. As the impurity potential is increased more 
and more induced charge is localized at these sites to more efficiently screen the 
impurity. 

For given magnitude of the impurity potential, we do not see a visual difference in 
the induced charge distribution between the positive and negative impurity. To be 
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Figure 3. Distribution of induced charges, on several lattice sites near the impurity, as a function of 
energy uo. Results are shown for increasing the spatial extent of the impurity potential by changing its 
Gaussian width cr/ao = 1.0 (a), 1.2 (b), 1.3 (c), and 1.4 (d), where ao is the lattice constant. We fixed 
C/o = 2.0t. Increasing the impurity size localizes more charges close to the impurity. It also enhances the 
weight of the induced charge on the next-nearest-neighbor lattice sites pointed out in the figure by arrows. 

more accurate about this visual inference, we quantitatively compare the plasmonic 
excitations at all energies as well as their corresponding spatial distribution of in- 
duced charges for positive and negative impurities (see Sec. 5 for more information 
on how we extract this information). Indeed, we confirm that our calculations sat- 
isfy the sum rule for all plasmonic excitations and that the corresponding induced 
charge distributions are identical for attractive and repelling impurity potential at 
half filling. 

The exact same response of the graphene lattice for both positive and negative 
single impurity at zero doping is related to the fundamental particle-hole symmetry 
of graphene. The tight-binding calculation shows that the valence and conduction 
band meet at two nonequivalent points in the Brillouin zone, which are called Dirac 
points. At zero doping the valence band is completely filled and the conduction 
band is empty (half filling). The Fermi level lies at the Dirac point so valence and 
conductance bands are symmetric with respect to the Fermi level (particle-hole 
symmetric). The same response of the lattice for positive and negative impurity 
potential originates from the particle-hole symmetric electronic band structure. 
We use this observation as a validity check for the numerical method implemented 
to calculate the plasmonic response in graphene. Later on we will see that once 
we break this symmetry, for example, by doping the system away from half filling 
(/X 7^ 0), the response of the system will be different for positive and negative single 
impurity. 

4.3. The effect of impurity size 

A spatially extended impurity is modeled by taking into account the effects of the 
impurity on lattice sites beyond nearest neighbor carbon sites. This is mathemat- 
ically achieved by relaxing the decay of the impurity potential over distance by 
increasing the magnitude of a/ao from 1.0 to 1.4, where ag is the lattice constant. 
Here we fix the potential to — 2.0t. 

In Figures [3]^ a)- (d), we show the effects of increased impurity size on the induced 
charge distributions. Comparing Fig. [sj^a) with Figs. [3]^b)-(d), one can see that as 
the size of the impurity potential is increased the localization of the induced charges 
in close vicinity of the impurity is increased (notice the increase in the range of the 
induced charge in the color legend of Fig. |3] compared to that of Fig. [2]). Moreover, 



August 23, 2011 0:8 Philosophical Magazine manuscript 'philmag'eps 



9 




Figure 4. Distribution of induced charges in graphene with positive impurity potential Uq = S.Ot and 
electron doping fi/t = 0.0 (a), 0.75 (b), 1.5 (c), and 2.5 (d). As the electron doping is increased the 
localization of the charges at sites close to the impurity site are enhanced, while it is decreased for lattice 
sites farther away. The number of localized plasmons (nano-plasmons) is increased as seen in the increasing 
number of bigger 3D sphere glyphs (charges) along the energy axis. In other words, screening by conduction 
electrons is increased in the vicinity of the single impurity. 

we see that it also enhances the weight of the induced charges on the next-nearest- 
neighbor carbon sites (compared to the magnitude of the charges at the lattice site 
where the arrow points to). Here too, we find that more and more induced charges 
are localized closer to the impurity with an increase of the extent of the impurity. 

4.4. The effect of positive and negative impurity potential at finite doping 

We model the effect of doping by taking a nonzero value for the chemical potential 
/i in Eq. ([T]). For electron doping the chemical potential moves above the Dirac 
point and lies in the conduction band, i.e., the chemical potential is positive. In 
Figures [4]^ a)- (d) we show the effects of the impurity on the induced charge distri- 
butions for doping levels /i/t = 0.0, 0.75, 1.5, and 2.5, respectively. The strength 
of the impurity potential is Uq = 3t. Comparing Fig. Qa) with Figs. |4]^b)-(d), 
one can see that as the level of electron doping is increased the localization of 
the induced charges at sites close to the impurity are enhanced. The magnitude 
of the induced charges on the remaining lattice sites is consequently suppressed. 
Moreover, the number of localized plasmon (nano-plasmons) increases as seen in 
the increased number of bigger 3D sphere glyphs along the energy axis. This effect 
is most clearly seen in Fig. [Ij 

The origin for the increase in the number of nano-plasmons (in the presence of a 
positively charged impurity), due to electron doping, comes from a change in the 
electronic and impurity states of graphene [21] . The impurity states lie at the band 
edge and electron doping brings the top of the filled graphene electronic states 
closer to the impurity states. This leads to an increase in the number of localized 
excitations near the impurity site. 

When the system is hole doped, the chemical potential moves below the Dirac 
point and lies in the valence band. The effect is modeled by taking a negative value 
for the chemical potential in Eq. Q. In Figure [5]^a)-(d) we show the effects of the 
impurity on the charge distributions for doping levels /j./t = 0.0, —0.75, —1.5, and 
—2.5, respectively. The strength of the impurity potential is fixed at Uq = 3.0t. 
Comparing Fig. [5]^ a) with the Figs. [5]^b)-(d), one can see that as the hole doping is 
increased the localization of charges in close vicinity of the impurity is suppressed. 
The suppression of localization of nano-plasmons is caused by the increased sepa- 
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Figure 5. Distribution of induced charges in graphene with positive impurity potential Uq = S.Ot and 
hole doping levels of /i/f = 0.0 (a), —0.75 (b), —1.5 (c), and —2.5 (d). As the hole doping is increased the 
localization of charges at sites close to the impurity is reduced. 




Figure 6. Probability function of induced charges over the entire 2D lattice in (a) pristine graphene 
and (b) graphene with a single impurity with a positive potential Uq = l.Ot. The distribution function of 
charges is symmetric with respect to the origin. The range of the magnitude of induced charges is increased 
in impure graphene. 



ration between the top of the fihed states of graphene and the impurity states. The 
different nano-plasmonic responses of the graphene lattice due to electron and hole 
doping originates from the breaking of particle hole-symmetry with doping away 
from half filling. 

To summarize, we studied nano-plasmonic excitations with a fixed negative im- 
purity potential but changing the chemical potential. We found that in the hole 
(electron) doped case the negative impurity potential leads to an increase (de- 
crease) of the localization of induced charges near the impurity. Thus a positive 
impurity potential for electron doping behaves similarly as a negative impurity po- 
tential for hole doping. Both of them enhance the localization around the impurity. 
In contrast, a positive impurity potential for hole doping and a negative impurity 
potential for electron doping suppress localization of induced charges around the 
impurity. These results are listed in Table [l] It is worth to notice that for a given 
magnitude of the chemical potential and impurity potential, flipping the sign of 
both of them gives identical nano-plasmonic response in energy and real space. 
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Figure 7. Induced charges in (a) pristine and (b) impure graphene as a function of energy at a fixed 
nearest-neighbor lattice site from the impurity, marked by the arrow in Fig. Qa). The result is for a 
positive impurity potential Uq = 2t. We see a larger range of induced charges in impure graphene. 



Table 1. Dependence of nano-plasmons in graphene on the sign of the impurity potential and doping (chemical 
potential). 
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Figure 8. The spatial distribution of the induced charges on every lattice site at given energy for pristine 
graphene. Only a few selected modes are shown. The energy w = uj/t oi plasmon labels the panels. 



5. Quantitative analysis 

Using an advanced visualization tool like ParaView one can extract in addition 
to the obvious visual cues [30j also quantitative information from the data set 
of interest. We demonstrate this capability by first plotting a histogram of the 
distribution of induced charges over all lattice sites and over all energies. The 
histogram is generated by applying the built-in 'Histogram' filter in ParaView [26] 
[27j on a data set with impurity potential Uq = l.Ot. This information may appear 
trivial, but is of interest to verify the conservation of total charge in our calculations, 
i.e., the sum of all induced charges must vanish and must be equally distributed 
between positive and negative oscillations. Since pristine graphene does not have 
localized plasmons, whereas impure graphene does, it is important to understand 
the rearrangement of induced charges. We can check for charge conservation and 
induced charge balance for two completely different cases, namely, localization vs. 
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Figure 9. The spatial distribution of the induced charges on individual lattice sites at given energy near 
the impurity on graphene lattice. Only a few selected nano-plasmons are shown. The energy w = uu/t of 
plasmon labels the panels. The localization of plasmons occurs at the nano-meter scale, hence the term 
nano-plasmons. 



delocalization of modes. The symmetric distribution of the induced charges as seen 
in the histogram or probabihty function of finding a specific induced charge in Fig.|6] 
imphes that the total induced charge is zero. We verified the charge neutrahty by 
calculating the sum of the induced charges. 

In both cases, pristine and impure graphene, the number of zero-induced charges 
is significantly larger than the number of nonzero induced charges, i.e, lattice sites 
that are unaffected. However the number of zero-induced charges is slightly bigger 
in impure graphene than in pristine graphene. The magnitude of induced charges 
decrease exponentially in impure graphene, while it has a Gaussian distribution 
in pristine graphene. The range of the magnitude of induced charges is slightly 
bigger for impure graphene. This implies that a relatively small number of plasmons 
acquires a bigger magnitude of induced charges. The larger number of zero-induced 
charges, as well as a wider range in the distribution of charges, and the exponential 
decrease are consequences of localized plasmonic excitations. 

Next, we monitor the induced charge at a single site in the graphene lattice as 
a function of energy. The calculation is performed with or without the presence of 
impurity. For the impurity case the impurity potential is fixed at Uq — 2.0t. We 
plot this information by using the built-in 'Plot Over Line' filter [261 ET]. It gives 
interpolated values of the variable along an arbitrary user-specified line. The result 
is shown in Fig. [Tlfor a particular nearest-neighbor site of the impurity, marked by 
the arrow in Fig. [l]^a). One also sees that the range of induced charges are bigger 
for the impure graphene. In addition, comparing Fig. |6]^b) and Fig. [7]^b) one can 
see that the range of induced charges in impure graphene increases with increasing 
impurity potential from — l.Ot to 2.0t. 

Finally, we extracted the energy-resolved spatial distribution of induced charges 
for pristine and impure graphene. As mentioned above, there are about 2300 plas- 
monic modes for a lattice with = 96 sites. It is very time consuming to extract 
every individual plasmonic excitation manually. However, this can be automated in 
ParaView by using the SUI (script user interface). We developed a simple python 
script, which controls the rendering in ParaView, to extract the distribution of 
the induced charges for all energies. It is worth mentioning that using a single- 
processor computer it is impractical to extract the necessary information in real 
time (it takes several hours to resolve the spatial distribution for all energies). This 
computationally intensive problem can be overcome by running the SUI in parallel 
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rendering mode. Since there are many unique plasmonic modes (characterized by 
the unique induced charge distributions), it is impossible to report ah of them. 
So we show only a few selected plasmonic modes in Fig. |8j One can clearly see 
that the cancellation of induced charges, which is required for the conservation of 
the total charge, can happen locally (e.g., pairwise) or globally (e.g., in groups). 
Some plasmons have simple patterns of dipolar, quadrupolar or triangular distri- 
bution of localized induced charges. However, others have non-trivial symmetries 
and correspond to higher-order multipoles. 

In Fig. [9] we show the energy-resolved spatial distribution of induced charges in 
impure graphene for a few selected plasmon energies with characteristic patterns. 
These charge distributions are very different from those of the pristine lattice. 
We can clearly see that the cancellation of induced charges can happen only over a 
larger area, because of the creation of localized plasmons. Some of these modes have 
simple patterns of dipolar, quadrupolar, radial or triangular symmetry. However, 
others have highly non-trivial patterns, e.g., clumps of similar charges on nearest 
lattice sites. The main result reported here is that the localization of plasmonic 
modes occurs at the nano-meter scale concentrating the electric field onto a very 
small area. 



6. Conclusions 

We investigated through exact diagonalization the plasmonic excitations of 
graphene in the presence of a single impurity with a special focus on localized 
modes (nano-plasmons) near the impurity. The dependence of excitations on the 
model parameters such as the sign, magnitude, size of the impurity potential and 
number of electrons in graphene was discussed. It was shown that the impurity 
potential and doping can be used to tune the properties of nano-plasmonic excita- 
tions, demonstrating that graphene is an inherently plasmonic material. 

Taking nano-plasmons in graphene as a case study, we demonstrated the impor- 
tance of a three-dimensional, interactive visualization tool to explore and analyze 
large data sets in a more effective and time efficient way. We described the visu- 
alization problem faced with, even when the data set is only several Megabytes in 
size. This issue arises when the data set is intrinsically multidimensional in both 
real space and parameter space. As a solution to this problem, we proposed a small- 
scale parallel processing approach, because it can no longer be efficiently handled 
by a single-processor computer. Finally, using an interactive visualization tool we 
could extract even more quantitative information leading to a better understanding 
of the physical properties of nano-plasmons on a two-dimensional lattice. 

The results discussed can be tested, in principle, in experiments by high- 
frequency optical probes or local scanning tunneling spectroscopy probes. However, 
in order to achieve such tests, calculations of extended or nano-sized impurity clus- 
ters are very desirable as they would give plasmonic excitations spread over a larger 
size with smaller frequencies. 
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